Accompanying code for annotating single-cell maps

Before beginning this tutorial, make sure you have the following packages installed:

R version	Packages
4.0.3	SingleCellExperiment_1.12.0*	dplyr_1.0.4*	ggplot2_3.3.3
.	Seurat_3.2.2	scmap_1.12.0*	Harmony_1.0**
.	scater_1.18.3*	celldex_1.0.0*	cerebroApp_1.3.0**
.	SCINA_1.2.0	SingleR_1.4.0*	msigdb_0.2.0**
.	devtools_2.3.2

"*" = packages must be installed by running BiocManager::install("package") instead of install.packages("package").
"**" = packages must be installed by from github using devtools, see code block below.

Alternatively, an installation code chunk has been provided below to help install the latest versions of all of the above packages. Run this only if you need to install all programs from scratch. Select compile from source if prompted.

This tutorial goes through multiple methods on annotation an unlabeled single-cell dataset, referred to as the query. We have chosen to label a single-cell experiment that consists of peripheral blood mononuclear cells (PBMCs). The tutorial consists of the following:

1. Reference-based automatic annotation This section annotates the query dataset using a previously labeled reference dataset. Many tools exist to do this: we are going over scmap (cell and cluster) and SingleR. We will further explore how to use integration as a form of annotation using Harmony.

2. Refining / Consensus annotations After finding multiple cell type labels for each cell using reference-based automatic annotation, we will keep the labels that most commonly occur across methods.

3. Marker-based automatic annotation Instead of using a reference dataset to annotate the query dataset, we will input lists of marker genes associated with specific cell types. The program we have chosen to demonstrate here is SCINA.

4. Manual annotation Here, we extract marker genes and associated pathways from the query dataset. To determine cell-type labels from this information, we would have to compare our differentially expressed genes and pathways to those described in the literature. To facilitate this process, we use Seurat and cerebroApp.

1. Reference-based automatic annotation

Create the Reference

The first step in performing reference-based annotation is to select an annotated dataset to use as the reference. Here we will use one of the references created by the authors of SingleR and show how it can be used with other tools such as scmap.

Other reference datasets can be found in GEO (https://www.ncbi.nlm.nih.gov/geo/) or at a link provided by the authors of the reference dataset. However, to use a dataset as a reference you will need both the single-cell RNA sequencing data and the cell-type annotations. GEO does not require authors to provide the cell-type annotations of their data, so you may need to contact the authors directly to to get the annotations for some datasets.

# Set a random seed to ensure result reproducibility
set.seed(9742)
# Download singleR reference data for immune cells and save it as the variable "ref"
# The variable is a class called "Summarized Experiment"
# This will take a while
ref <- celldex::DatabaseImmuneCellExpressionData()

Next we need to reformat the data to ensure it is compatible with the tool we are using. We will be demonstrating scmap, which uses data formatted as a ‘SingleCellExperiment object’, and assumes by default that gene names are found in a column named ‘feature_symbol’ while the cell-type labels are in a column named ‘cell_type1’. In addition, scmap requires that you normalize and log-transform the reference data; this has already been done for the SingleR reference data so we skip those steps here.

# Assign cell-type labels in a column named "cell_type1"
colData(ref)$cell_type1 <- colData(ref)$label.fine
# Assign gene names in a column called "feature_symbol"
rowData(ref)$feature_symbol <- rownames(ref)

# Convert the data into a SingleCellExperiment object
ref_sce <- SingleCellExperiment::SingleCellExperiment(assays=list(logcounts=Matrix::Matrix(assays(ref)$logcounts)), 
            colData=colData(ref), rowData=rowData(ref))

Our reference data is ready to be used now. So lets process this data to build the index we will use to map our unlabeled data to. First, we select genes to use, which will be those deemed most informative by scmap after fitting a linear model to the gene expression by gene dropout distribution. Those which are most informative have high expression values and low % dropout rates across cells.

# Create scmap-cluster reference by first selecting the most informative features
ref_sce <- scmap::selectFeatures(ref_sce, suppress_plot=FALSE)

Your object does not contain counts() slot. Dropouts were calculated using logcounts() slot...

# Inspect the first 50 genes selected by scmap
rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)][1:50]

 [1] "ABCD2"      "AC006129.2" "AC008372.1" "AC009831.1" "AC010760.1" "AC013264.2" "AC018867.1" "AC067945.4"
 [9] "AC068641.1" "AC079150.2" "AC087239.1" "AC092580.4" "AC093110.3" "AC093642.4" "AC103563.9" "ACVR1C"    
[17] "ADTRP"      "AL109767.1" "AL356776.1" "AL450992.2" "ALOX5AP"    "AMICA1"     "AMIGO1"     "ANXA1"     
[25] "AP000304.1" "APBA2"      "AQP3"       "ATP8B2"     "AXIN2"      "BACH2"      "BCAS4"      "BCL11B"    
[33] "BEX2"       "BEX5"       "BTLA"       "C10orf35"   "C14orf64"   "C20orf197"  "CACNA1I"    "CAMK4"     
[41] "CARD11"     "CCDC64"     "CCL5"       "CCR4"       "CCR7"       "CD2"        "CD200"      "CD226"     
[49] "CD247"      "CD27"      

# You can check and see how many genes were chosen by checking the length of the
# vector of gene names
length(rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)])

[1] 500

Now we can see the genes that scmap has chosen to use. If there are key marker genes missing we can make sure they are included like this:

# Create a list of key markers that you want to use
my_key_markers = c("TRAC", "TRBC1", "TRBC2", "TRDC", "TRGC1", "TRGC2", "IGKC")
# Ensure markers are in the list of features used by scmap
rowData(ref_sce)$scmap_features[rownames(ref_sce) %in% my_key_markers] <- TRUE
# You can check and see if this added any genes by checking the length 
# of the vector of gene names again
length(rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)])

[1] 502

And we can remove genes that we think might be technical artefacts, such as mitochondria RNAs, like this:

# Create a list of mitochondrial genes from the dataset (genes that begin with "MT")
mt_genes <- rownames(ref_sce)[grep("^MT-", rownames(ref_sce))]
# Remove these genes from the features used by scmap
rowData(ref_sce)$scmap_features[rownames(ref_sce) %in% mt_genes] <- FALSE
# Check how many genes this is
length(rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)])

[1] 495

# Extract the features and assign them to a new variable, "scmap_feature_genes"
scmap_feature_genes <- rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)]
# Note that the number of genes/features is identical to what we just checked
length(scmap_feature_genes)

[1] 495

Now we build the reference profiles used in scmap-cluster, for cluster-based cell-type annotation. These profiles can be accessed and plotted from inside the SingleCellExperiment object as follows:

# Create reference profiles;
# Once reference profiles are generated the original data are 
# not needed for scmap-cluster
ref_sce <- scmap::indexCluster(ref_sce)
# Visualize interesting features as a heatmap
# Reformat the data so that they can be used as input to ggplot2
cormat <- reshape2::melt(as.matrix(metadata(ref_sce)$scmap_cluster_index))
# Plot the data
ggplot2::ggplot(cormat, ggplot2::aes(x = Var2, y = Var1, fill = value)) +
  ggplot2::geom_tile() +
  ggplot2::scale_fill_gradient2(low = "blue", high = "darkred",
                                name = "Expression value") +
  ggplot2::theme_minimal() +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 1,
                                            size = 18, hjust = 1),
                 axis.text.y = ggplot2::element_text(size = 15),
                 axis.title.x = ggplot2::element_blank(),
                 axis.title.y = ggplot2::element_blank())

# Store expression information as a variable
scmap_cluster_reference <- metadata(ref_sce)$scmap_cluster_index

From here on out, scmap only needs this set of reference profiles. So if working with a very large reference, one could save this index separately to your computer and reload it when annotating new datasets. But since that is not the case here, we will simply save this index to a variable for now.

We will also demonstrate scmap-cell to annotate individual cells of our dataset, so we will create that index as well. As before one would first normalize and log-transform the reference data, and select genes to use. As we have already done that, we need only run the command to build the scmap-cell index. There are two parameters we can set: M and k, increasing M and k will give more accurate mapping but increase the size of the index, and the time needed to map cells. Here we use the defaults (you may see a warning message about the defaults that are being used):

# Update the previous reference to also contain the scmap-cell reference
ref_sce <- scmap::indexCell(ref_sce)

Parameter M was not provided, will use M = n_features / 10 (if n_features <= 1000), where n_features is the number of selected features, and M = 100 otherwise.
Parameter k was not provided, will use k = sqrt(number_of_cells)

# Extract the scmap index from the reference and store as a variable
scmap_cell_reference <- metadata(ref_sce)$scmap_cell_index
# Extract the associated cell IDs from the reference and save as a variable
scmap_cell_metadata <- colData(ref_sce)

scmap-cell assigns cells in one dataset to their “nearest neighbours” in the reference dataset. In this case, the “nearest neighbours” are the cells in the reference dataset most similar to the cells in the query dataset.

One can use any rule they like to transfer information, such as cell-type or pseudotime, from these nearest neighbours to the query data. Thus we need to store the associated metadata (cell type ID) for the reference as well (see above). Now we don’t need to use our original reference dataset anymore.

Assign cells from the query dataset to the reference.

The query dataset we will be using is provided by 10X genomics.

download.file("https://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz",
              "pbmc3k_filtered_gene_bc_matrices.tar.gz")

trying URL 'https://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz'
Content type 'application/x-tar' length 7621991 bytes (7.3 MB)
==================================================
downloaded 7.3 MB

untar("pbmc3k_filtered_gene_bc_matrices.tar.gz")

Now we need to load our unlabeled dataset into R. Normal preprocessing including QC filtering, normalizing and log-transforming the data must be done prior to annotating. In addition, scmap is based on the SingleCellExperiment object, so if our data is stored as a Seurat object we must convert it to SingleCellExperiment as shown below.

# This portion of the tutorial is assuming the raw 10X data is in the
# following folder in your directory:
data <- Seurat::Read10X("filtered_gene_bc_matrices/hg19/")
# Make SingleCellExperiment from the raw matrix
query_sce <- SingleCellExperiment::SingleCellExperiment(assays=list(counts=data))

# Make SingleCellExperiment from Seurat object
query_seur <- Seurat::CreateSeuratObject(data)

Feature names cannot have underscores ('_'), replacing with dashes ('-')

query_sce <- Seurat::as.SingleCellExperiment(query_seur)

# normalize the data using the scater package
query_sce <- scater::logNormCounts(query_sce)

# add feature_symbol column (i.e. the gene symbols)
rowData(query_sce)$feature_symbol <- rownames(query_sce)

Now you should have an entry in assays(my_sce) called logcounts with the log-normalized matrix. We are now ready to annotate our data with scmap-cluster. Let’s start with scmap-cluster:

# Run scmapCluster
scmap_cluster_res <- scmap::scmapCluster(projection=query_sce, 
                index_list=list(immune1 = scmap_cluster_reference), 
                threshold=0.1)

Features AC009831.1, AC067945.4, AC068641.1, AC079150.2, AL356776.1, CHRM3-AS2, CLUHP3, CTB-118N6.1, IGHV1OR15-1, IGKC, LA16c-390H2.4, MINOS1P3, MIR199A1, MIR221, MIR3136, MIR3176, MIR320E, MIR3671, MIR378I, MIR4257, MIR4292, MIR4420, MIR4489, MIR4520A, MIR4639, MIR4645, MIR4697, MIR4712, MIR4742, MIR548AR, MIR5581, MIR5685, MIR641, PSPHP1, RNA5SP37, RNA5SP82, RNU1-106P, RNU1-60P, RNU5A-6P, RNU6-322P, RNU6-415P, RNU6-5P, RNU6-878P, RNU6ATAC18P, RNU7-140P, RNU7-163P, RNU7-49P, RNU7-57P, RNU7-75P, RP11-104L21.2, RP11-1094M14.4, RP11-1094M14.5, RP11-1094M14.8, RP11-1094M14.9, RP11-23J18.1, RP11-317B3.2, RP11-332E4.1, RP11-402J7.2, RP11-436I9.6, RP11-564A8.8, RP11-603B24.1, RP11-686D22.10, RP11-832A4.5, RP4-575N6.5, RP4-631H13.6, RP4-765C7.2, RPL12P18, RPL32P1, RPL39P, RPS4XP22, RPSAP53, SCARNA18, SNORA14, SNORA19, SNORA26, SNORD101, SNORD104, SNORD108, SNORD12B, SNORD14E, SNORD42B, SNORD51, SNORD6, SNORD60, SNORD72, SNORD78, SNORD82, SNORD99, SOD1P3, TRAC, TRAJ1, TRAJ10, TRAJ11, TRAJ12, TRAJ13, TRAJ14, TRAJ15, TRAJ16, TRAJ17, TRAJ18, TRAJ19, TRAJ2, TRAJ20, TRAJ21, TRAJ22, TRAJ23, TRAJ24, TRAJ25, TRAJ26, TRAJ27, TRAJ28, TRAJ29, TRAJ3, TRAJ30, TRAJ31, TRAJ32, TRAJ33, TRAJ34, TRAJ35, TRAJ36, TRAJ37, TRAJ38, TRAJ39, TRAJ4, TRAJ40, TRAJ41, TRAJ42, TRAJ43, TRAJ44, TRAJ45, TRAJ47, TRAJ48, TRAJ49, TRAJ5, TRAJ52, TRAJ53, TRAJ54, TRAJ57, TRAJ6, TRAJ7, TRAJ8, TRAJ9, TRAV1-1, TRAV1-2, TRAV10, TRAV12-1, TRAV12-2, TRAV12-3, TRAV13-1, TRAV13-2, TRAV14DV4, TRAV16, TRAV17, TRAV18, TRAV19, TRAV2, TRAV20, TRAV21, TRAV22, TRAV23DV6, TRAV24, TRAV25, TRAV26-1, TRAV26-2, TRAV27, TRAV29DV5, TRAV3, TRAV30, TRAV35, TRAV36DV7, TRAV38-1, TRAV38-2DV8, TRAV39, TRAV4, TRAV41, TRAV5, TRAV6, TRAV8-1, TRAV8-2, TRAV8-3, TRAV8-4, TRAV8-6, TRAV9-2, TRBC2, TRBJ2-2, TRBJ2-3, TRBJ2-7, TRBV1, TRBV10-2, TRBV11-1, TRBV19, TRBV2, TRBV20-1, TRBV21-1, TRBV23-1, TRBV24-1, TRBV25-1, TRBV27, TRBV28, TRBV29-1, TRBV3-1, TRBV30, TRBV4-1, TRBV4-2, TRBV5-1, TRBV5-4, TRBV5-5, TRBV5-6, TRBV6-1, TRBV6-5, TRBV6-6, TRBV7-3, TRBV7-6, TRBV7-8, TRBV9, TRDC, TRDV1, TRGC1, TRGC2, TRGJP2, TRGV10, TRGV3, TRGV4, Z95152.1 are not present in the 'SCESet' object and therefore were not set.

# plot the results of our annotation
par(mar=c(13, 4, 1, 0))
barplot(table(scmap_cluster_res$combined_labs), las=2)

# Store this annotation information within the query object
colData(query_sce)$scmap_cluster <- scmap_cluster_res$combined_labs

# Make a UMAP of the cells, labeled with the cell-type annotations from scmapCluster
query_sce <- scater::runUMAP(query_sce)
scater::plotReducedDim(query_sce, dimred="UMAP", colour_by="scmap_cluster")

Alternatively we could use scmap-cell, to find the 10 nearest neighbours to each cell (i.e. the 10 most similar cells to each query cell), then pick the annotation that is most common among the neighbours, like this:

# Determine the 10 nearest neighbours from the reference dataset for each
# cell in the query dataset using scmapCell
nearest_neighbours <- scmap::scmapCell(projection=query_sce, 
    index_list = list(immune1 = scmap_cell_reference), 
    w=10)

# Get metadata (cell type IDs) for the neighbours of each cell in the query dataset
mode_label <- function(neighbours, metadata=scmap_cell_metadata$cell_type1) {
    freq <- table(metadata[neighbours])
    label <- names(freq)[which(freq == max(freq))]
    if (length(label) > 1) {return("ambiguous")}
    return(label)
}

# Apply these labels to the query cells
scmap_cell_labs <- apply(nearest_neighbours$immune1$cells, 2, mode_label)

# Add the labels to the query object
colData(query_sce)$scmap_cell <- scmap_cell_labs

# Create a bar plot of how many cells in the query dataset were assigned
# a specific label
par(mar=c(10, 4, 0, 0))
barplot(table(scmap_cell_labs), las=2)

# Make a UMAP and add the new cell-type annotations
scater::plotReducedDim(query_sce, dimred="UMAP", colour_by="scmap_cell")

Another option compatible with the SingleCellExperiment Object is SingleR. As before, we need a reference and a query dataset. In the case of SingleR, we need the entirety of the reference dataset, rather than generating a compressed reference index as we did with scmap. In addition, running just this small example demonstrates the difference in run time between the methods (SingleR takes a fair bit of time).

# Run SingleR on the query data and the reference to acquire
# cell-type predictions for the cells in the query dataset
predictions <- SingleR::SingleR(test=query_sce, ref=ref, labels=ref$label.fine)

Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 106Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 106Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 105Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 105Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 104Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 103Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 104Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 102Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 104Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 104Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 104Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 104Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 104Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 104Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 102Argument 'center' should be of the same length as number of rows of 'x'. Use of a scalar value is deprecated: 1 != 2700

# You'll notice that some of the cells didn't get assigned a cell identity
# We can count the number here:
sum(is.na(predictions$pruned.labels))

[1] 23

# Change NAs to "ambiguous"
predictions$pruned.labels[which(is.na(predictions$pruned.labels))] <- "ambiguous"
# Add singleR labels to query_sce
colData(query_sce)$singleR <- predictions$pruned.labels

# Create a bar plot of number of cells per assigned cell ID
par(mar=c(13, 4, 2, 0))
barplot(table(predictions$pruned.labels), las=2)

# Make a UMAP and add the cell-type annotations
scater::plotReducedDim(query_sce, dimred="UMAP", colour_by="singleR")

Integration as a form of annotation

Another option is to integrate our query data with our reference data. Then we simply transfer the labels from the annotated reference to the neighbouring query cells in the integrated dataset. Clustering the integrated data is a common approach to transferring labels. We demonstrate how this could be done with Harmony below. But the approach would be the same for any integration tool.

Note: the SingleR reference is not single cells, but averages across many cells. Thus we convert and downsample the reference to a single cell object for demonstration purposes. For a real experiment, one would use the original single cells as the reference when integrating datasets.

set.seed(2891)
# Convert reference and query datasets to Seurat Objects

# Add a "counts" slot to the reference SingleCellExperiment object so we can convert it to a Seurat Object
assays(ref_sce)[["counts"]] <- round(2^assays(ref_sce)[["logcounts"]]) -1
colnames(ref_sce) <- paste("cell", 1:ncol(ref_sce))

# Subset both objects so both the reference and query datasets have the same genes
# First subset the reference
ref_seur <- Seurat::as.Seurat(ref_sce[rownames(ref_sce) %in% rownames(query_sce),])
ref_seur@active.ident <- factor(rep("reference", ncol(ref_seur)))
# Now subset the query
query_seur <- Seurat::as.Seurat(query_sce[rownames(query_seur) %in% rownames(ref_sce),])
query_seur@active.ident <- factor(rep("query", ncol(query_seur)))

# Downsample the reference to be similar to query in terms of total UMIs
totalUMI <- median(query_seur@meta.data$nCount_RNA)
ref_seur@assays$RNA@counts <- Seurat::SampleUMI(ref_seur@assays$RNA@counts,
                                                max.umi=totalUMI, upsample=FALSE)

# Merge the datasets together into a single Seurat object
merged_seur <- merge(ref_seur, query_seur)
merged_seur@meta.data$source <- merged_seur@active.ident

# Normalize the combined data
merged_seur <- Seurat::NormalizeData(merged_seur)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

# Rather than choosing new variable features, we will choose
# the genes that had been previously important by scmap for consistency
Seurat::VariableFeatures(merged_seur) <- scmap_feature_genes

Not all features provided are in this Assay object, removing the following feature(s): AC009831.1, AC067945.4, AC068641.1, AC079150.2, AL356776.1, CHRM3-AS2, CLUHP3, CTB-118N6.1, IGHV1OR15-1, IGKC, LA16c-390H2.4, MINOS1P3, MIR199A1, MIR221, MIR3136, MIR3176, MIR320E, MIR3671, MIR378I, MIR4257, MIR4292, MIR4420, MIR4489, MIR4520A, MIR4639, MIR4645, MIR4697, MIR4712, MIR4742, MIR548AR, MIR5581, MIR5685, MIR641, PSPHP1, RNA5SP37, RNA5SP82, RNU1-106P, RNU1-60P, RNU5A-6P, RNU6-322P, RNU6-415P, RNU6-5P, RNU6-878P, RNU6ATAC18P, RNU7-140P, RNU7-163P, RNU7-49P, RNU7-57P, RNU7-75P, RP11-104L21.2, RP11-1094M14.4, RP11-1094M14.5, RP11-1094M14.8, RP11-1094M14.9, RP11-23J18.1, RP11-317B3.2, RP11-332E4.1, RP11-402J7.2, RP11-436I9.6, RP11-564A8.8, RP11-603B24.1, RP11-686D22.10, RP11-832A4.5, RP4-575N6.5, RP4-631H13.6, RP4-765C7.2, RPL12P18, RPL32P1, RPL39P, RPS4XP22, RPSAP53, SCARNA18, SNORA14, SNORA19, SNORA26, SNORD101, SNORD104, SNORD108, SNORD12B, SNORD14E, SNORD42B, SNORD51, SNORD6, SNORD60, SNORD72, SNORD78, SNORD82, SNORD99, SOD1P3, TRAC, TRAJ1, TRAJ10, TRAJ11, TRAJ12, TRAJ13, TRAJ14, TRAJ15, TRAJ16, TRAJ17, TRAJ18, TRAJ19, TRAJ2, TRAJ20, TRAJ21, TRAJ22, TRAJ23, TRAJ24, TRAJ25, TRAJ26, TRAJ27, TRAJ28, TRAJ29, TRAJ3, TRAJ30, TRAJ31, TRAJ32, TRAJ33, TRAJ34, TRAJ35, TRAJ36, TRAJ37, TRAJ38, TRAJ39, TRAJ4, TRAJ40, TRAJ41, TRAJ42, TRAJ43, TRAJ44, TRAJ45, TRAJ47, TRAJ48, TRAJ49, TRAJ5, TRAJ52, TRAJ53, TRAJ54, TRAJ57, TRAJ6, TRAJ7, TRAJ8, TRAJ9, TRAV1-1, TRAV1-2, TRAV10, TRAV12-1, TRAV12-2, TRAV12-3, TRAV13-1, TRAV13-2, TRAV14DV4, TRAV16, TRAV17, TRAV18, TRAV19, TRAV2, TRAV20, TRAV21, TRAV22, TRAV23DV6, TRAV24, TRAV25, TRAV26-1, TRAV26-2, TRAV27, TRAV29DV5, TRAV3, TRAV30, TRAV35, TRAV36DV7, TRAV38-1, TRAV38-2DV8, TRAV39, TRAV4, TRAV41, TRAV5, TRAV6, TRAV8-1, TRAV8-2, TRAV8-3, TRAV8-4, TRAV8-6, TRAV9-2, TRBC2, TRBJ2-2, TRBJ2-3, TRBJ2-7, TRBV1, TRBV10-2, TRBV11-1, TRBV19, TRBV2, TRBV20-1, TRBV21-1, TRBV23-1, TRBV24-1, TRBV25-1, TRBV27, TRBV28, TRBV29-1, TRBV3-1, TRBV30, TRBV4-1, TRBV4-2, TRBV5-1, TRBV5-4, TRBV5-5, TRBV5-6, TRBV6-1, TRBV6-5, TRBV6-6, TRBV7-3, TRBV7-6, TRBV7-8, TRBV9, TRDC, TRDV1, TRGC1, TRGC2, TRGJP2, TRGV10, TRGV3, TRGV4, Z95152.1

# Scale the data and run dimensionality reduction on the combined data
merged_seur <- Seurat::ScaleData(merged_seur)

Centering and scaling data matrix

  |                                                                                                         
  |                                                                                                   |   0%
  |                                                                                                         
  |===================================================================================================| 100%

merged_seur <- Seurat::RunPCA(merged_seur)

PC_ 1 
Positive:  CD28, RCAN3, ICOS, PKIA, CD96, IL7R, CD2, TRIB2, CTLA4, ITK 
       BCL11B, TC2N, SLAMF1, PPP1R16B, ATP8B2, IL2RA, GPRIN3, GIMAP6, LEF1, HAPLN3 
       FAM101B, CD3D, SH2D1A, RP11-326C3.2, HPCAL4, CD40LG, CD3E, IL6R, CCR7, C14orf64 
Negative:  LYZ, CST3, FCER1G, FCN1, S100A9, S100A11, S100A8, IFITM3, CYBB, CD300E 
       CSTA, SYK, MT2A, FCGR3A, RPS4Y1, FAM49A, SMCO4, S100Z, IFNGR2, ANXA1 
       IL8, SECTM1, AMICA1, TNFRSF1A, GIMAP7, NCF4, SEPT1, GZMM, PRF1, HRH2 
PC_ 2 
Positive:  IL32, LCK, CD3E, GZMA, CCL5, KLRB1, CST7, CD247, GZMM, TC2N 
       SKAP1, TRIB2, MAL, PRF1, TIGIT, PYHIN1, AQP3, RPS4Y1, SLA2, GPRIN3 
       ATP8B2, CD3D, BCL11B, ALOX5AP, GZMK, CD27, SEPT1, LINC00861, PRKCQ-AS1, C14orf64 
Negative:  TNFSF14, TNFRSF9, CD200, IL2, PSAT1, NDFIP2, LIF, BTLA, IL2RA, BACH2 
       CCDC64, IL8, IER3, IFNG, PKIA, SLAMF1, RGS16, YES1, GPR171, IFNGR2 
       CD40LG, ICOS, TBX21, FCN1, CYBB, S100A9, RP11-536K7.5, CSF1, CD300E, CSTA 
PC_ 3 
Positive:  CD300E, CSTA, GNAQ, FCN1, CYBB, S100A9, CD4, FAM101B, CST3, FCER1G 
       LYZ, IFITM3, FAM49A, S100A8, S100Z, SYK, S100A11, IL6R, TRIB2, SMCO4 
       IFNGR2, GPRIN3, IL8, SECTM1, ATP8B2, LY96, RCAN3, SULT1B1, CD28, NCF4 
Negative:  IL32, CD7, GIMAP7, CD3D, TNFSF14, TNFRSF9, RPS4Y1, GZMM, CD200, GIMAP5 
       PSAT1, NDFIP2, IL2, RHOH, IFNG, CD3E, TBX21, GPR171, OCIAD2, CCL5 
       LIF, CD27, CD8A, SEPT1, GZMK, GZMA, CD3G, CST7, BTLA, CCDC64 
PC_ 4 
Positive:  LEF1, CCR7, CD3D, MAL, RPS4Y1, C14orf64, IL7R, SUSD3, CD27, RCAN3 
       TRAT1, NGFRAP1, CAMK4, AQP3, RHOH, PRKCQ-AS1, ADTRP, CD3E, EPHX2, TSHZ2 
       TCEA3, CD3G, EDAR, CD28, IL6R, PRKCA, INPP4B, LINC00402, GNG7, OCIAD2 
Negative:  CST7, GZMA, PRF1, CCL5, IL2RB, TBX21, FCGR3A, TGFBR3, KLRB1, CD247 
       PYHIN1, SH2D2A, CD226, FCER1G, IFNG, GZMM, SH2D1A, SBK1, S100B, PLEKHG3 
       IL18R1, C20orf197, CTD-2562J17.7, CD8A, FAM49A, SLA2, RP11-473M20.9, PPP1R16B, TSPAN2, ALOX5AP 
PC_ 5 
Positive:  ANXA1, S100A11, IL32, DUSP4, CD4, CTLA4, CCR4, GIMAP7, GIMAP4, CD3D 
       CD40LG, AQP3, LYZ, MT2A, TIGIT, S100A9, GATA3, S100A8, CST3, POLR2L 
       FCN1, CD3E, GPR171, CD2, IL2RA, IL6R, TRAT1, AMICA1, ICOS, NPDC1 
Negative:  GNG7, BACH2, IL24, CD8A, S100B, KRT72, SYK, SH3BGRL2, BTLA, NELL2 
       CCR7, RRAS2, IKZF3, EDAR, C14orf64, SULT1B1, BCAS4, PLXDC1, FAM229B, FAM101B 
       HOOK1, CYBB, P2RY10, AC018867.1, AC093110.3, AC087239.1, AMIGO1, THEMIS, ZNF793, RP11-75L1.1 

merged_seur <- Seurat::RunUMAP(merged_seur, dims=1:15)

The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session07:33:57 UMAP embedding parameters a = 0.9922 b = 1.112
07:33:57 Read 4261 rows and found 15 numeric columns
07:33:57 Using Annoy for neighbor search, n_neighbors = 30
07:33:57 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
07:33:58 Writing NN index file to temp file /var/folders/n9/vs4199cs1dd63v34z3gqn1dc0000gn/T//Rtmp8VECRS/fileca3312236a8a
07:33:58 Searching Annoy index using 1 thread, search_k = 3000
07:33:59 Annoy recall = 100%
07:34:00 Commencing smooth kNN distance calibration using 1 thread
07:34:02 Found 2 connected components, falling back to 'spca' initialization with init_sdev = 1
07:34:02 Initializing from PCA
07:34:02 PCA: 2 components explained 47.6% variance
07:34:02 Commencing optimization for 500 epochs, with 179004 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
07:34:08 Optimization finished

Seurat::DimPlot(merged_seur, reduction="umap") + ggplot2::ggtitle("Before Integration")

# Run Harmony to remove batch effects
merged_seur <- harmony::RunHarmony(merged_seur, "source", dims.use=1:15)

Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 8/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 9/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 10/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.RNA.harmony; see ?make.names for more details on syntax validity

merged_seur <- Seurat::RunUMAP(merged_seur, dims=1:15, reduction="harmony")

07:34:13 UMAP embedding parameters a = 0.9922 b = 1.112
07:34:13 Read 4261 rows and found 15 numeric columns
07:34:13 Using Annoy for neighbor search, n_neighbors = 30
07:34:13 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
07:34:13 Writing NN index file to temp file /var/folders/n9/vs4199cs1dd63v34z3gqn1dc0000gn/T//Rtmp8VECRS/fileca33223d9135
07:34:13 Searching Annoy index using 1 thread, search_k = 3000
07:34:14 Annoy recall = 100%
07:34:15 Commencing smooth kNN distance calibration using 1 thread
07:34:17 Found 3 connected components, falling back to 'spca' initialization with init_sdev = 1
07:34:17 Initializing from PCA
07:34:17 PCA: 2 components explained 53.34% variance
07:34:17 Commencing optimization for 500 epochs, with 183288 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
07:34:23 Optimization finished

# Plot the data
Seurat::DimPlot(merged_seur, reduction="umap") + ggplot2::ggtitle("After Integration")

Now that the data is integrated we will cluster the data and look at the annotations of the reference cells present in each cluster. As with all clustering, this may require manual tuning of the resolution parameters to get the best labels.

# Cluster the integrated dataset
merged_seur <- Seurat::FindNeighbors(merged_seur, reduction="harmony", dims=1:15)

Computing nearest neighbor graph
Computing SNN

merged_seur <- Seurat::FindClusters(merged_seur, resolution=0.5)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 4261
Number of edges: 187562

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9200
Number of communities: 10
Elapsed time: 0 seconds

# Plot the data
Seurat::DimPlot(merged_seur, reduction="umap") + ggplot2::ggtitle("After Integration")

# Create a table of cluster labels based on integrated data
table(merged_seur@meta.data$label.fine, 
        merged_seur@active.ident)

                                  
                                     0   1   2   3   4   5   6   7   8   9
  B cells, naive                     0   0   0 106   0   0   0   0   0   0
  Monocytes, CD14+                   0   0 106   0   0   0   0   0   0   0
  Monocytes, CD16+                   0   0   0   0   0   0 105   0   0   0
  NK cells                           0 105   0   0   0   0   0   0   0   0
  T cells, CD4+, memory TREG         0   0   0   0   0   0   0 104   0   0
  T cells, CD4+, naive               4   0   0   0  95   4   0   0   0   0
  T cells, CD4+, naive TREG          1   0   0   0 102   0   0   1   0   0
  T cells, CD4+, naive, stimulated   0   0   0   0   0   0   0   0 102   0
  T cells, CD4+, TFH                96   0   0   0   5   0   0   3   0   0
  T cells, CD4+, Th1               103   0   0   0   1   0   0   0   0   0
  T cells, CD4+, Th1_17            104   0   0   0   0   0   0   0   0   0
  T cells, CD4+, Th17               99   0   0   0   0   0   0   5   0   0
  T cells, CD4+, Th2                96   0   0   0   6   0   0   2   0   0
  T cells, CD8+, naive               0   0   0   0   1 103   0   0   0   0
  T cells, CD8+, naive, stimulated   0   0   0   0   0   1   0   0   0 101

Here we have a table of the reference annotations (across rows) per cluster (across columns). We can manually label the clusters based on this table or we could create a rule to algorithmically label the clusters based on this table. Since there are only 11 clusters, we assign the labels manually.

cluster_labs <- c("0"="ambiguous", 
    "1"="Monocytes, CD14+", 
    "2"="B cells, naive", 
    "3"="T cells, CD4+, naive TREG",
    "4"="T cells, CD4+, Th1_17",
    "5"="NK cells",
    "6"="T cells, CD8+, naive",
    "7"="Monocytes, CD16+",
    "8"="T cells, CD4+, memory TREG",
    "9"="T cells, CD4+, naive, stimulated",
    "10" = "T cells, CD8+, naive, stimulated")

# Assign cluster label to the associated query cells
# (the query cells that had been assigned the same cluster label)
merged_seur@meta.data$annotation <- cluster_labs[merged_seur@meta.data$RNA_snn_res.0.5]

# Add the results to the SingleCellExperiment Object and plot
query_sce$Harmony_lab <- merged_seur@meta.data$annotation[merged_seur@meta.data$source =="query"]
scater::plotReducedDim(query_sce, dimred="UMAP", colour_by="Harmony_lab")

2. Refining / Consensus annotations

Once we have run several tools, we can use the consensus of the labels to get a more robust annotation. In this case we will simply use the most common label across tools to assign the final automatically annotated label.

annotation_columns <- c("scmap_cluster", "scmap_cell", "singleR", "Harmony_lab")

#Optional check how consistent the labelling was.
#head(colData(query_sce)[,annotation_columns])

get_consensus_label <- function(labels){
    labels <- labels[labels != "ambiguous"]
    if (length(labels) == 0) {return("ambiguous")}
    freq <- table(labels)
    label <- names(freq)[which(freq == max(freq))]
    if (length(label) > 1) {return("ambiguous")}
    return(label)
}

colData(query_sce)$consensus_lab <- apply(colData(query_sce)[,annotation_columns], 1, get_consensus_label)
scater::plotReducedDim(query_sce, dimred="UMAP", colour_by="consensus_lab")

3. Marker-based automatic annotation

An alternative way for annotation of your query scRNAseq dataset is to utilize Marker-based annotation tools. SCINA is a semi-supervised annotation tool that takes in the signature genes and expression matrix and predicts the potential labels based on the prior knowledge of the cell-type-specific markers. List of markers is usually provided in the gmt format. The PBMC gene set used below have been gathered by Diaz-Mejia JJ et al.

download.file("https://zenodo.org/record/3369934/files/pbmc_22_10x.tar.bz2",
              "pbmc_22_10x.tar.bz2")

trying URL 'https://zenodo.org/record/3369934/files/pbmc_22_10x.tar.bz2'
Content type 'application/octet-stream' length 374541 bytes (365 KB)
==================================================
downloaded 365 KB

untar("pbmc_22_10x.tar.bz2")

The extracted data will by located in the following file: ./MY_PAPER/SUPPLEMENTARY_DATA/pbmc_22_10x/pbmc_22_10x_cell_type_signature_gene_sets.gmt

The results from this annotation tool are not used in the above step to find consensus annotations because the lists of marker genes are not consistent with the cell types identified in the reference dataset. This is because these data come from different sources, and would not have been characterizing the exact same set of cells. If you wish for marker-based and reference-based annotation methods to be combined in the above step of automatically determining consensus annotations, you would have to make sure all of the identified cell subtypes are the same and that they are spelt the exact same way in order for R to recognize the names as identical.

# Import the marker genes as a GMT file and store as a variable
markers <- msigdb::read.gmt('./MY_PAPER/SUPPLEMENTARY_DATA/pbmc_22_10x/pbmc_22_10x_cell_type_signature_gene_sets.gmt')
# Convert the expression data from Seurat object into a matrix data structure
exprMatrix <- as.matrix(Seurat::GetAssayData(query_seur))
# Run SCINA on the query data using the marker genes to identify cell types
# Specifying rm_overlap = FALSE allows the same marker gene to specify multiple cell types which
# may be useful if identifying cell subtypes or other similar types of cells
# Specifying allow_unknown = TRUE allows cells to be labeled as "unknown" instead of being
# assigned a low-confident label
predictions.scina = SCINA::SCINA(exp = exprMatrix, signatures = markers$genesets,
                          rm_overlap = FALSE, allow_unknown = TRUE)
# Add SCINA annotation information to each cell in Seurat object
colData(query_sce)$SCINA <- predictions.scina$cell_labels

# Make a UMAP and add the SCINA cell-type annotations
scater::plotReducedDim(query_sce, dimred="UMAP", colour_by="SCINA") +
  ggplot2::theme(legend.position = "bottom",
                 legend.text = ggplot2::element_text(size = 4))

4. Manual annotation

Retrieving marker genes

If you do not have an extensive list of markers per cell type, or a good quality reference dataset, it is useful to extract the top marker genes from each cluster of your query data. We can easily do this in Seurat, with the data formatted as a *Seurat object (which we created earlier and stored as the variable query_seur). First, the data must be normalized and scaled, and the variable genes between cells must be determined.

query_seur <- Seurat::NormalizeData(query_seur) # Normalize the data

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

query_seur <- Seurat::FindVariableFeatures(query_seur) # Determine the variable features of the dataset

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

query_seur <- Seurat::ScaleData(query_seur) # Scale the data based on the variable features

Centering and scaling data matrix

  |                                                                                                         
  |                                                                                                   |   0%
  |                                                                                                         
  |==================================================                                                 |  50%
  |                                                                                                         
  |===================================================================================================| 100%

Next, different types of dimensionality reduction must be performed on the data so that the cells can be grouped together in 2D space.

query_seur <- Seurat::RunPCA(query_seur)

PC_ 1 
Positive:  CST3, TYROBP, LST1, AIF1, FTL, FCN1, LYZ, FTH1, S100A9, FCER1G 
       TYMP, CFD, LGALS1, S100A8, LGALS2, CTSS, SERPINA1, SPI1, IFITM3, PSAP 
       CFP, SAT1, IFI30, COTL1, S100A11, NPC2, LGALS3, GSTP1, PYCARD, NCF2 
Negative:  MALAT1, LTB, IL32, CD2, ACAP1, STK17A, CTSW, CD247, CCL5, GIMAP5 
       AQP3, GZMA, CST7, TRAF3IP3, MAL, HOPX, ITM2A, GZMK, MYC, GIMAP7 
       BEX2, ETS1, LDLRAP1, ZAP70, LYAR, RIC3, TNFAIP8, KLRG1, NKG7, SAMD3 
PC_ 2 
Positive:  CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DRA, HLA-DQB1, LINC00926, CD79B, HLA-DRB1, CD74 
       HLA-DPB1, HLA-DMA, HLA-DQA2, HLA-DRB5, HLA-DPA1, HLA-DMB, FCRLA, HVCN1, LTB, BLNK 
       KIAA0125, P2RX5, IRF8, IGLL5, SWAP70, ARHGAP24, SMIM14, PPP1R14A, FCRL2, C16orf74 
Negative:  NKG7, PRF1, CST7, GZMA, GZMB, FGFBP2, CTSW, GNLY, GZMH, SPON2 
       CCL4, FCGR3A, CCL5, CD247, XCL2, CLIC3, AKR1C3, SRGN, HOPX, CTSC 
       TTC38, S100A4, ANXA1, IL32, IGFBP7, ID2, ACTB, XCL1, APOBEC3G, SAMD3 
PC_ 3 
Positive:  HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1, CD74, HLA-DPA1, MS4A1, HLA-DRB1, HLA-DRB5 
       HLA-DRA, HLA-DQA2, TCL1A, LINC00926, HLA-DMB, HLA-DMA, HVCN1, FCRLA, IRF8, BLNK 
       KIAA0125, SMIM14, PLD4, P2RX5, IGLL5, SWAP70, TMSB10, LAT2, IGJ, MZB1 
Negative:  PPBP, PF4, SDPR, SPARC, GNG11, NRGN, GP9, RGS18, TUBB1, CLU 
       HIST1H2AC, AP001189.4, ITGA2B, CD9, TMEM40, CA2, PTCRA, ACRBP, MMD, NGFRAP1 
       TREML1, F13A1, RUFY1, SEPT5, MPP1, TSC22D1, CMTM5, RP11-367G6.3, MYL9, GP1BA 
PC_ 4 
Positive:  VIM, S100A8, S100A6, S100A4, TMSB10, S100A9, IL32, GIMAP7, S100A10, LGALS2 
       RBP7, MAL, FCN1, LYZ, CD2, S100A12, MS4A6A, FYB, S100A11, AQP3 
       GIMAP4, FOLR3, ANXA1, MALAT1, AIF1, GIMAP5, IL8, IFI6, TRABD2A, ASGR1 
Negative:  HLA-DQA1, CD79A, CD79B, HIST1H2AC, PF4, SDPR, HLA-DQB1, PPBP, MS4A1, CD74 
       GNG11, SPARC, HLA-DPB1, HLA-DQA2, GP9, HLA-DRB1, HLA-DPA1, NRGN, RGS18, HLA-DRA 
       TCL1A, PTCRA, CD9, AP001189.4, CA2, CLU, TUBB1, LINC00926, HLA-DRB5, ITGA2B 
PC_ 5 
Positive:  GZMB, FGFBP2, NKG7, GNLY, PRF1, CCL4, CST7, SPON2, GZMA, CLIC3 
       GZMH, XCL2, CTSW, TTC38, CCL5, AKR1C3, IGFBP7, XCL1, S100A8, CCL3 
       TYROBP, HOPX, CD160, HAVCR2, S100A9, PTGDR, FCER1G, LGALS2, RBP7, S100A12 
Negative:  LTB, VIM, AQP3, PPA1, MAL, KIAA0101, CD2, CYTIP, CORO1B, FYB 
       IL32, TRADD, ANXA5, HN1, TUBA1B, TYMS, PTGES3, ITM2A, COTL1, GPR183 
       ACTG1, TNFAIP8, ATP5C1, TRAF3IP3, ZWINT, GIMAP4, PRDX1, ABRACL, LDLRAP1, NGFRAP1 

query_seur <- Seurat::RunTSNE(query_seur)
# RunUMAP has already been performed on the data, so the following line of code
# does not need to be run in this case:
#query_seur <- Seurat::RunUMAP(query_seur, dims = 1:50)

From this object, we can cluster the data at a chosen resolution that can be modified later on if desired.

# Determine the "nearest neighbours" of each cell
query_seur <- Seurat::FindNeighbors(query_seur, dims = 1:50)

Computing nearest neighbor graph
Computing SNN

# Cluster the cells
query_seur <- Seurat::FindClusters(query_seur, resolution = 0.5)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 2700
Number of edges: 205239

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8313
Number of communities: 7
Elapsed time: 0 seconds

Before extracting the marker genes, let’s visualize our data on a UMAP.

Seurat::DimPlot(query_seur, reduction = "UMAP")

Now let’s extract the top marker genes, and see which ones correspond with each cluster. This can be done using the FindAllMarkers function within Seurat.

markers_seur <- Seurat::FindAllMarkers(query_seur, only.pos = TRUE)

# Markers are now stored and can be viewed in the following table. They are ordered from lowest to highest p-value percluster:
markers_seur

The expression of marker genes across clusters are commonly viewed as dot plots or heat maps. The dot plot communicates the percentage of cells in a cluster a marker gene is expressed in (the size of the dot) and mean detected gene expression for that cluster. A heat map communicates average marker gene expression across clusters. Both graphs are created below using the top 5 marker genes detected by Seurat for each cluster (filtered by the package dplyr) and default settings for the plots.

require(dplyr)

Loading required package: dplyr

Attaching package: ‘dplyr’

The following object is masked from ‘package:Biobase’:

    combine

The following objects are masked from ‘package:GenomicRanges’:

    intersect, setdiff, union

The following object is masked from ‘package:GenomeInfoDb’:

    intersect

The following objects are masked from ‘package:IRanges’:

    collapse, desc, intersect, setdiff, slice, union

The following object is masked from ‘package:matrixStats’:

    count

The following objects are masked from ‘package:S4Vectors’:

    first, intersect, rename, setdiff, setequal, union

The following objects are masked from ‘package:BiocGenerics’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

# Retrieve the top 5 marker genes per cluster
# Use whichever genes have the highest values under the AVG_LOG column
top5 <- markers_seur %>% group_by(cluster) %>%
  dplyr::slice_max(get(grep("^avg_log", colnames(markers_seur), value = TRUE)),
                   n = 5)
# Create the dot plot
Seurat::DotPlot(query_seur, features = unique(top5$gene)) +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 1,
                                            size = 8, hjust = 1)) +
  Seurat::NoLegend()

# Create the heatmap
Seurat::DoHeatmap(query_seur, features = unique(top5$gene)) +
  Seurat::NoLegend() +
  ggplot2::theme(axis.text.y = ggplot2::element_text(size = 8))

The following features were omitted as they were not found in the scale.data slot for the RNA assay: NOSIP, CD3E, CD3D, IL7R, LDHB

Pathway analysis

Pathway analysis can also be done for each cluster to determine significantly up- and downregulated pathways based on known gene function. An easy way to do this is by feeding our current Seurat object into cerebroApp. cerebroApp requires that marker genes be fetched again through before performing simple pathway analysis.

# First get marker genes through cerebro
query_seur <- cerebroApp::getMarkerGenes(query_seur,
                                         groups = c('seurat_clusters'),
                                         assay = "RNA",
                                         organism = "hg")

# Get enriched pathways through cerebro
query_seur <- cerebroApp::getEnrichedPathways(query_seur,
                                              databases = c("GO_Biological_Process_2018",
                                                            "GO_Cellular_Component_2018",
                                                            "GO_Molecular_Function_2018",
                                                            "KEGG_2016",
                                                            "WikiPathways_2016",
                                                            "Reactome_2016",
                                                            "Panther_2016",
                                                            "Human_Gene_Atlas",
                                                            "Mouse_Gene_Atlas"),
                                              adj_p_cutoff = 0.05,
                                              max_terms = 100,
                                              URL_API = "http://amp.pharm.mssm.edu/Enrichr/enrich")

# Enriched pathways are stored in the following location:
query_seur@misc$enriched_pathways

$cerebro_seurat_enrichr
$cerebro_seurat_enrichr$seurat_clusters

NANA

Combining top marker genes with functional pathway information should be fairly indicative of cell type depending on your selected clustering resolution. You can easily rerun the analyses after modifying the resolution to get a better idea of subtypes (increase resolution) or identify more general cell types (decrease resolution).

Pathway analysis can also be done by pasting a list of marker genes for a specific cluster into an online resource such as gProfiler.

None of these methods can determine the identity of a cell with absolute certainty. However, combining these resources can provide robust support for cell-type labels in a query dataset, allowing for a confidently labeled single-cell map.